# Compute profit

	function compute_plan_profit(plan_total_premiums,average_claims,plan_transfers,plan_demands,var_profit_flag)
		claims_reins = (1 .- pmt[:,findall(pmt_fields .== "reins_factor")[1]]) .* average_claims .* plan_demands;
		variable_admin_costs = pmt[:,findall(pmt_fields .== "variable_admin")[1]] .* plan_demands
		plan_profit = transpose(foc_moment_indicator) * (plan_total_premiums .- claims_reins .+ 
			plan_transfers .- variable_admin_costs) ./
			(transpose(foc_moment_indicator) * plan_demands) .- average_plan_fixed_costs  .- foc_residuals; 
		return(plan_profit)
	end

# Tabulate coverage

	function tabulate_coverage_by_metal(plan_demands)
		bronze_coverage = sum(plan_demands[plans_pmt[!,:AV] .== .6])/total_weight * total_pop;
		silver_coverage = sum(plan_demands[plans_pmt[!,:AV] .== .7])/total_weight * total_pop;
		gold_coverage = sum(plan_demands[plans_pmt[!,:AV] .== .8])/total_weight * total_pop;
		platinum_coverage = sum(plan_demands[plans_pmt[!,:AV] .== .9])/total_weight * total_pop;
		total_coverage = bronze_coverage + silver_coverage + gold_coverage + platinum_coverage;
		return(cat(dims=1,bronze_coverage,silver_coverage,gold_coverage,platinum_coverage,total_coverage))
	end
	
	function tabulate_coverage_by_network_type(plan_demands)
		HMO_coverage = sum(plan_demands[plans_pmt[!,:HMO] .== 1])/total_weight * total_pop;
		PPO_coverage = sum(plan_demands[plans_pmt[!,:HMO] .== 0])/total_weight * total_pop;
		return(cat(dims=1,HMO_coverage,PPO_coverage))
	end
	
	function tabulate_coverage_by_insurer(plan_demands)
		coverage = zeros(length(insurers));
		for f in 1:length(insurers)
			coverage[f] = sum(plan_demands[plans_pmt[!,:Insurer] .== insurers[f]])/total_weight * total_pop;
		end
		return(coverage)
	end
	
	function get_dominated_enrollment(probability_weights)
		gold_dominated = sum(probability_weights .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .8) .* (households[household_numbers,:FPL] .<= 2));
		platinum_dominated = sum(probability_weights .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .9) .* (households[household_numbers,:FPL] .<= 1.5));
		return(cat(dims=1,0,gold_dominated + platinum_dominated,-gold_dominated,-platinum_dominated))
	end
	
# Compute Premiums	
	
	function compute_avg_premium_by_metal(plan_total_premiums,plan_demands)
		
		avs = [.6,.7,.8,.9];
		average_premiums = zeros(length(avs));
		for m in 1:length(avs)
			metal_indices = findall(plans_pmt[!,:AV] .== avs[m]);
			average_premiums[m] = sum(plan_total_premiums[metal_indices])/sum(plan_demands[metal_indices])
		end
		
		return(average_premiums)
	end
	
	function compute_avg_premium_by_firm(plan_total_premiums,plan_demands)
		
		average_premiums = zeros(length(insurers));
		for f in 1:length(insurers)
			insurer_indices = findall(plans_pmt[!,:Insurer] .== insurers[f]);
			average_premiums[f] = sum(plan_total_premiums[insurer_indices])/sum(plan_demands[insurer_indices])
		end
		
		return(average_premiums)
	end
	
	function compute_avg_premium_by_nettype(plan_total_premiums,plan_demands)
	
		avg_HMO_premium = sum(plan_total_premiums[findall(plans_pmt[!,:HMO] .== 1)])/
			sum(plan_demands[findall(plans_pmt[!,:HMO] .== 1)]);
		avg_PPO_premium = sum(plan_total_premiums[findall(plans_pmt[!,:HMO] .== 0)])/
			sum(plan_demands[findall(plans_pmt[!,:HMO] .== 0)]);
		
		return(cat(dims=1,avg_HMO_premium,avg_PPO_premium))
	end
	
# Compute Claims	
	
	function compute_avg_claims_by_metal(average_claims,average_pt_claims,plan_demands)
	
		avs = [.6,.7,.8,.9];
		average_claims_by_metal = zeros(length(avs));
		for m in 1:length(avs)
			metal_indices = findall(plans_pmt[!,:AV] .== avs[m]);
			average_claims_by_metal[m] = sum(plan_demands[metal_indices] .* average_claims[metal_indices]) /
				sum(plan_demands[metal_indices]);
		end
		
		return(average_claims_by_metal)
	end
	
	function compute_avg_claims_by_insurer(average_claims,average_pt_claims,plan_demands)
		average_claims_by_insurer = zeros(length(insurers));
		for f in 1:length(insurers)
			average_claims_by_insurer[f] = sum(plan_demands[plans_pmt[!,:Insurer] .== insurers[f]] .* 
				average_claims[plans_pmt[!,:Insurer] .== insurers[f]])/sum(plan_demands[plans_pmt[!,:Insurer] .== insurers[f]]);
		end
		return(average_claims_by_insurer)
	end
	
# FOC residuals	
	
	function calculate_foc_resids_by_metal(plan_pt_demands,foc_resids)
	
		avs = [.6,.7,.8,.9];
		foc_resids_by_metal = zeros(length(avs));
		for m in 1:length(avs)
			metal_indices = findall(plans[!,:AV] .== avs[m]);
			foc_resids_by_metal[m] = 
				sum(plan_pt_demands[metal_indices] .* foc_resids[metal_indices])/sum(plan_pt_demands[metal_indices])
		end
		
		return(foc_resids_by_metal)
	end
	
	function calculate_foc_resids_by_firm(plan_pt_demands,foc_resids)
	
		foc_resids_by_firm = zeros(length(insurers));
		for f in 1:length(insurers)
			insurer_indices = findall(plans[!,:Insurer] .== insurers[f]);
			foc_resids_by_firm[f] = 
				sum(plan_pt_demands[insurer_indices] .* foc_resids[insurer_indices])/sum(plan_pt_demands[insurer_indices])
		end
	
		return(foc_resids_by_firm)
	end
	
# Output Processing


	output_fields = ["Consumer Surplus";"Profit";"Premium Subsidies";
						"Cost Sharing Subsidies";"Penalties";"Uncompensated Care";"Social Welfare";
						"Consumer Surplus - Kappa 75%";"Consumer Surplus - Kappa 50%";"Consumer Surplus - Kappa 25%";"Consumer Surplus - Kappa 0%";
						"Bronze";"Silver";"Gold";"Platinum";"Total Coverage";"% Enrolled";
						"% Switching Plan (Base)";"% Switching Metal (Base)";"% Switching Insurer (Base)";
						"% Switching Plan (Setting)";"% Switching Metal (Setting)";"% Switching Insurer (Setting)";
						"HMO";"PPO";
						insurers;
						"Bronze Premium";"Silver Premium";"Gold Premium";"Platinum Premium";
						"Anthem Premium";"Blue Shield Premium";"Health Net Premium";"Kaiser Premium";"Other Insurer Premium";
						"HMO Premium";"PPO Premium";"Average Premium";"Average Subsidized Premium";
						"Bronze Premium (Intensive Selection)";"Silver Premium (Intensive Selection)";
						"Gold Premium (Intensive Selection)";"Platinum Premium (Intensive Selection)";"Average Premium (Intensive Selection)";
						"Bronze Premium (No Metal Selection)";"Silver Premium (No Metal Selection)";
						"Gold Premium (No Metal Selection)";"Platinum Premium (No Metal Selection)";"Average Premium (No Metal Selection)";
						"Bronze Premium (No Selection)";"Silver Premium (No Selection)";
						"Gold Premium (No Selection)";"Platinum Premium (No Selection)";"Average Premium (No Selection)";
						"Bronze Avg Claims Resid";"Silver Avg Claims Resid";"Gold Avg Claims Resid";"Platinum Avg Claims Resid";"Average Claims Resids";
						"Bronze Avg Claims (Intensive Selection)";"Silver Avg Claims (Intensive Selection)";"Gold Avg Claims (Intensive Selection)";"Platinum Avg Claims (Intensive Selection)";"Average Claims (Intensive Selection)";
						"Bronze Avg Claims (No Metal Selection)";"Silver Avg Claims (No Metal Selection)";"Gold Avg Claims (No Metal Selection)";"Platinum Avg Claims (No Metal Selection)";"Average Claims (No Metal Selection)";
						"Bronze Avg Claims (No Selection)";"Silver Avg Claims (No Selection)";"Gold Avg Claims (No Selection)";"Platinum Avg Claims (No Selection)";"Average Claims (No Selection)";
						"Anthem Avg Claims Resid";"Blue Shield Avg Claims Resid";"Health Net Avg Claims Resid";"Kaiser Avg Claims Resid";"Other Insurer Avg Claims Resid";
						"Anthem Avg Claims (No Selection)";"Blue Shield Avg Claims (No Selection)";"Health Net Avg Claims (No Selection)";"Kaiser Avg Claims (No Selection)";"Other Insurer Avg Claims (No Selection)";
						"Uninsured Risk Score";"Average Risk Score";"Average AV";"Average FOC Residuals";"Max FOC Residuals";
						"CWlt250";"CW250to400";"CWgt400";"CW0to17";"CW18to34";"CW35to54";"CW55plus";
						"CWFemale";"CWMale";"CWSingle";"CWFamily";"CWAsian";"CWBlack";"CWHispanic";"CWWhite";"CWOther"];
		
	output = zeros(length(output_fields));	
		
	total_weight = sum(households[!,:weight])
	
	# Total population is sum of household weights (which are 5% sample) and flagged household weight
	# which was made by make.data.objects.R
	cd(datadir)
	flagged_market_size_by_year = CSV.read("flagged_market_size_by_year.csv",DataFrame,header=true); 
	total_pop = total_weight/.05 + flagged_market_size_by_year[flagged_market_size_by_year[!,:year] .== year,:market_size][1];
		
	# Get Base probabilities
	
		if choice_seq_flag
			if scenario == "ACA"
				base_base_premium = base_premium;
			elseif network_flag == true
				base_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			else 
				base_base_premium = CSV.read(string(year,"Base0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			end
		
		elseif all_rs_variables
			if scenario == "ACA"
				base_base_premium = base_premium;
			elseif network_flag == true
				base_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			else 
				base_base_premium = CSV.read(string(year,"Base0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			end
		elseif int_rs_variables
			if scenario == "ACA"
				base_base_premium = base_premium;
			elseif inattention_flag
				base_base_premium = CSV.read(string(year,"Inattention0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif network_flag == true
				base_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			else 
				base_base_premium = CSV.read(string(year,"Base0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			end
		else
			if scenario == "ACA"
				base_base_premium = base_premium;
			elseif network_flag == true
				base_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_noseq.csv"),DataFrame,header=true)[!,:x3];
			else 
				base_base_premium = CSV.read(string(year,"Base0_premium_output_noseq.csv"),DataFrame,header=true)[!,:x3];
			end
		end
		
		# Initialize key variables
		base_subsidy = copy(voucher);
		base_new_unsubs_premiums = base_base_premium[household_pt_indices] .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		X_base = copy(X);
		if inattention_flag
			Xsearch_base = copy(Xsearch);
		end
		
		# Create data frame to hold objects we will group by household
		base_silverdf = DataFrame(household_number = household_numbers[silver_indicators],
			premium = base_new_unsubs_premiums[silver_indicators],base_silver_indices = silver_indicators);
		
		# Update Subsidy
		
			# Determine SLC plans by household
			base_num_household_silver_plans = combine(groupby(base_silverdf,:household_number), :household_number => length)[!,:household_number_length];
			pop!(base_num_household_silver_plans)
			base_num_household_silver_plans = cat(dims=1,0,base_num_household_silver_plans);
			base_starting_indices = cumsum(base_num_household_silver_plans); # starting position for each household
			base_household_slc_indices = combine(groupby(base_silverdf,:household_number), :premium => determine_SLC_plan)[!,:premium_determine_SLC_plan];
			base_slc_indices = base_silverdf[base_starting_indices .+ base_household_slc_indices,:base_silver_indices]; # SLC indices in hdata
			base_slc_plans = zeros(size(hdata)[1]);
			base_slc_plans[base_slc_indices] = (households[!,:subsidized_members] .> 0); # SLC index gets 1 if house is subsidized, 0 otherwise
		
			# Determine new subsidy 
			base_subsidy = max.(0,base_new_unsubs_premiums[base_slc_indices] - households[!,:SLC_contribution]) .* 
				(households[!,:subsidized_members] .> 0); # each household's subsidy (0 if unsubsidized);
			base_subsidy = base_subsidy[household_numbers];
						
		# Update Premium Variables
		X_base[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]] = 
				(max.(0,base_new_unsubs_premiums .- base_subsidy) .- 
				(1 - mand_repeal_flag) .* hdata[:,findall(hdata_fields .== "penalty")[1]]) ./ hdata[:,findall(hdata_fields .== "members")[1]];			
		X_base[uninsured_plans .== 0,premium_indices] = 
				X_base[uninsured_plans .== 0,premium_indices] .* X_base[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]];
		if random_parameters
			X_base[:,random_parameter_indices] = X_base[:,random_parameter_indices] .* V[:,:,1]; 
		end
		
		if inattention_flag & (year > 2014)
		
			# Growth rates
			growth_rates = base_new_unsubs_premiums ./ forward_simulated_prev_premiums[uninsured_plans .== 0];
			fsim_prev_choices = forward_simulated_prev_choices[uninsured_plans .== 0];
			growthdf = DataFrame(household_number = household_numbers[offered_indices],growth_rate = growth_rates[offered_indices],prev_choice=fsim_prev_choices[offered_indices]);
			median_growth_rates = combine(groupby(growthdf,:household_number), :growth_rate => median);
			
			household_median_growth_rates = zeros(size(households)[1]);
			household_median_growth_rates[median_growth_rates[!,:household_number]] = median_growth_rates[!,:growth_rate_median];
			
			growthdf_for_choices = growthdf[findall(growthdf[!,:prev_choice] .== 1),:];
			growthdf_for_choices[!,:median_growth_rate] = household_median_growth_rates[growthdf_for_choices[!,:household_number]];
			
			premium_shock_households = growthdf_for_choices[findall(growthdf_for_choices[!,:growth_rate] .> growthdf_for_choices[!,:median_growth_rate]),:household_number];
			Xsearch_base[:,findall(search_covariates .== "premium_shock")[1]] .= 0;
			Xsearch_base[premium_shock_households,findall(search_covariates .== "premium_shock")[1]] .= 1;
		
		end
		
		base_expdf = DataFrame(household_number = household_numbers,
					expvarall = exp.((X_base[uninsured_plans .== 0,:] * min_beta)/lambda));
		#if except_network_flag
		#	base_expdf[!,:expvarall] = exp.((X_base[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
		#end

		# Update choice probabilities
		base_expsum = combine(groupby(base_expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
		if inattention_flag & (year > 2014)
			entrants = findall(households[!,:previous_plan_number] .== "NA");
			search_probabilities = 1 ./ (1 .+ exp.(Xsearch_base * beta_search));
			search_probabilities[entrants] .= 1; 
					
			base_uninsured_probabilities = 1  ./ (1 .+ base_expsum);
			base_probabilities = search_probabilities[household_numbers] .* base_expdf[!,:expvarall] ./ (1 .+ base_expsum[household_numbers]) .+ 
				(1 .- search_probabilities[household_numbers]) .* data[uninsured_plans .== 0,:previous_choice];
		else
			base_uninsured_probabilities = 1  ./ (1 .+ (base_expsum).^(lambda));
			base_probabilities = base_expdf[!,:expvarall] .* (base_expsum[household_numbers]).^(lambda-1) ./ (1 .+ (base_expsum[household_numbers]).^(lambda));
		end
		
		base_probability_weights = base_probabilities .* hdata[:,findall(hdata_fields .== "weight")[1]];
		
			# ideally I'd get rid of these ACA variables
			aca_probabilities = copy(base_probabilities); 
			aca_probability_weights = copy(base_probability_weights);
			aca_plan_demands = zeros(length(plan_pmt_names));
			for j in 1:length(plan_pmt_names)	
				hh_plan_indices = findall(household_pmt_indices .== j);
				pt_index = findall(plans[!,:pt_name] .== year_plan_names[j])[1];
				aca_plan_demands[j] = sum(aca_probability_weights[hh_plan_indices]);	
			end
		
####### NOTE: aca_average_Claims needs to be fixed to the Base case, not ACA
	# This affects the "No Sorting" average claims panel (avg claims should be the same across all scenarios in this panel)
		
			X_CLMS[:,findall(claims_moments_covariates .== "log_risk_score")[1]] = X_RS * theta_RS;
			aca_average_claims = exp.(X_CLMS * theta_CLMS);	
		
		base_all_probabilities = zeros(size(X_base)[1]);
		base_all_probabilities[uninsured_plans .== 0] = base_probabilities;
		base_all_probabilities[uninsured_plans .== 1] = base_uninsured_probabilities;
		base_all_probability_weights = base_all_probabilities .* data[!,:weight];
		
		base_im_weight_object = DataFrame(household_number = household_numbers_data,
			bronze_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== .6),
			silver_weights = base_all_probability_weights .* ((X_base[:,findall(all_covariates .== "AV")[1]] .== .7) .|
														   (X_base[:,findall(all_covariates .== "AV")[1]] .== .73) .|
														   (X_base[:,findall(all_covariates .== "AV")[1]] .== .87) .|
														   (X_base[:,findall(all_covariates .== "AV")[1]] .== .94)),
			gold_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== .8),
			platinum_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== .9),
			uninsured_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== 0),
			anthem_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "Anthem")[1]] .== 1),
			bs_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "Blue_Shield")[1]] .== 1),
			hn_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "Health_Net")[1]] .== 1),
			kaiser_weights = base_all_probability_weights .* (X_base[:,findall(all_covariates .== "Kaiser")[1]] .== 1),
			small_weights = base_all_probability_weights .* (1 .- (X_base[:,findall(all_covariates .== "Anthem")[1]] .+
				X_base[:,findall(all_covariates .== "Blue_Shield")[1]] .+
				X_base[:,findall(all_covariates .== "Health_Net")[1]] .+ 
				X_base[:,findall(all_covariates .== "Kaiser")[1]])) .*  
				(X_base[:,findall(all_covariates .== "AV")[1]] .> 0));
		
		base_im_household_probability_weights = 
			DataFrame(Bronze = combine(groupby(base_im_weight_object,:household_number), :bronze_weights => sum)[!,:bronze_weights_sum],
					  Silver = combine(groupby(base_im_weight_object,:household_number), :silver_weights => sum)[!,:silver_weights_sum],
					  Gold = combine(groupby(base_im_weight_object,:household_number), :gold_weights => sum)[!,:gold_weights_sum],
					  Platinum = combine(groupby(base_im_weight_object,:household_number), :platinum_weights => sum)[!,:platinum_weights_sum],
					  Uninsured = combine(groupby(base_im_weight_object,:household_number), :uninsured_weights => sum)[!,:uninsured_weights_sum],
					  Anthem = combine(groupby(base_im_weight_object,:household_number), :anthem_weights => sum)[!,:anthem_weights_sum],
					  Blue_Shield = combine(groupby(base_im_weight_object,:household_number), :bs_weights => sum)[!,:bs_weights_sum],
					  Health_Net = combine(groupby(base_im_weight_object,:household_number), :hn_weights => sum)[!,:hn_weights_sum],
					  Kaiser = combine(groupby(base_im_weight_object,:household_number), :kaiser_weights => sum)[!,:kaiser_weights_sum],
					  Small_Insurer = combine(groupby(base_im_weight_object,:household_number), :small_weights => sum)[!,:small_weights_sum]);
		
	# Get Base Setting probabilities (Same as whatever we are running, but without inertia)
		
		if choice_seq_flag
			if (pc_flag & nora_flag)
				setting_base_premium = CSV.read(string(year,"No_RA_PC0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif pc_flag
				setting_base_premium = CSV.read(string(year,"PC0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif nora_flag
				setting_base_premium = CSV.read(string(year,"No_RA0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif nochurn_flag
				setting_base_premium = CSV.read(string(year,"No_Churn0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif eq_robust_flag
				setting_base_premium = CSV.read(string(year,"Robust0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif voucher_flag
				setting_base_premium = CSV.read(string(year,"Voucher0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif network_flag
				setting_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif network_flag
				setting_base_premium = CSV.read(string(year,"Inattention0_premium_output_intvars_seq_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif scenario == "ACA"
				setting_base_premium = base_premium;
			else
				setting_base_premium = base_base_premium;
			end
		elseif all_rs_variables
			if (pc_flag & nora_flag)
				setting_base_premium = CSV.read(string(year,"No_RA_PC0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif pc_flag
				setting_base_premium = CSV.read(string(year,"PC0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif nora_flag
				setting_base_premium = CSV.read(string(year,"No_RA0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif nochurn_flag
				setting_base_premium = CSV.read(string(year,"No_Churn0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif eq_robust_flag
				setting_base_premium = CSV.read(string(year,"Robust0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif voucher_flag
				setting_base_premium = CSV.read(string(year,"Voucher0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif network_flag
				setting_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_allvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif scenario == "ACA"
				setting_base_premium = base_premium;
			else
				setting_base_premium = base_base_premium;
			end
		else
			if (pc_flag & nora_flag)
				setting_base_premium = CSV.read(string(year,"No_RA_PC0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif pc_flag
				setting_base_premium = CSV.read(string(year,"PC0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif nora_flag
				setting_base_premium = CSV.read(string(year,"No_RA0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif nochurn_flag
				setting_base_premium = CSV.read(string(year,"No_Churn0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif eq_robust_flag
				setting_base_premium = CSV.read(string(year,"Robust0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif voucher_flag
				setting_base_premium = CSV.read(string(year,"Voucher0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif network_flag
				setting_base_premium = CSV.read(string(year,"ACA_Net0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif inattention_flag
				setting_base_premium = CSV.read(string(year,"Inattention0_premium_output_intvars_avint.csv"),DataFrame,header=true)[!,:x3];
			elseif scenario == "ACA"
				setting_base_premium = base_premium;
			else
				setting_base_premium = base_base_premium;
			end
		end
	
	
		# Initialize key variables
		setting_subsidy = copy(voucher);
		setting_new_unsubs_premiums = setting_base_premium[household_pt_indices] .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		X_base = copy(X);
		if inattention_flag
			Xsearch_base = copy(Xsearch);
		end
		
		# Create data frame to hold objects we will group by household
		setting_silverdf = DataFrame(household_number = household_numbers[silver_indicators],premium = setting_new_unsubs_premiums[silver_indicators],
			setting_silver_indices = silver_indicators);
		
		# Update Subsidy
		
			# Determine SLC plans by household
			setting_num_household_silver_plans = combine(groupby(setting_silverdf,:household_number), :household_number => length)[!,:household_number_length];
			pop!(setting_num_household_silver_plans)
			setting_num_household_silver_plans = cat(dims=1,0,setting_num_household_silver_plans);
			setting_starting_indices = cumsum(setting_num_household_silver_plans); # starting position for each household
			setting_household_slc_indices = combine(groupby(setting_silverdf,:household_number), :premium => determine_SLC_plan)[!,:premium_determine_SLC_plan];
			setting_slc_indices = setting_silverdf[setting_starting_indices .+ setting_household_slc_indices,:setting_silver_indices]; # SLC indices in hdata
			setting_slc_plans = zeros(size(hdata)[1]);
			setting_slc_plans[setting_slc_indices] = (households[!,:subsidized_members] .> 0); # SLC index gets 1 if house is subsidized, 0 otherwise
		
			# Determine new subsidy 
			if !voucher_flag
				setting_subsidy = max.(0,setting_new_unsubs_premiums[setting_slc_indices] - households[!,:SLC_contribution]) .* 
					(households[!,:subsidized_members] .> 0); # each household's subsidy (0 if unsubsidized);
				setting_subsidy = setting_subsidy[household_numbers];
			end		

		# Update Premium Variables
		X_base[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]] = 
				(max.(0,setting_new_unsubs_premiums .- setting_subsidy) .- 
				(1 - mand_repeal_flag) .* hdata[:,findall(hdata_fields .== "penalty")[1]]) ./ hdata[:,findall(hdata_fields .== "members")[1]];			
		X_base[uninsured_plans .== 0,premium_indices] = 
				X_base[uninsured_plans .== 0,premium_indices] .* X_base[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]];
		if random_parameters
			X_base[:,random_parameter_indices] = X_base[:,random_parameter_indices] .* V[:,:,1]; 
		end
		
		if inattention_flag & (year > 2014)
		
			# Growth rates
			growth_rates = setting_new_unsubs_premiums ./ forward_simulated_prev_premiums[uninsured_plans .== 0];
			fsim_prev_choices = forward_simulated_prev_choices[uninsured_plans .== 0];
			growthdf = DataFrame(household_number = household_numbers[offered_indices],growth_rate = growth_rates[offered_indices],prev_choice=fsim_prev_choices[offered_indices]);
			median_growth_rates = combine(groupby(growthdf,:household_number), :growth_rate => median);
			
			household_median_growth_rates = zeros(size(households)[1]);
			household_median_growth_rates[median_growth_rates[!,:household_number]] = median_growth_rates[!,:growth_rate_median];
			
			growthdf_for_choices = growthdf[findall(growthdf[!,:prev_choice] .== 1),:];
			growthdf_for_choices[!,:median_growth_rate] = household_median_growth_rates[growthdf_for_choices[!,:household_number]];
			
			premium_shock_households = growthdf_for_choices[findall(growthdf_for_choices[!,:growth_rate] .> growthdf_for_choices[!,:median_growth_rate]),:household_number];
			Xsearch_base[:,findall(search_covariates .== "premium_shock")[1]] .= 0;
			Xsearch_base[premium_shock_households,findall(search_covariates .== "premium_shock")[1]] .= 1;
		
		end
		
		setting_expdf = DataFrame(household_number = household_numbers,
					expvarall = exp.((X_base[uninsured_plans .== 0,:] * min_beta)/lambda));
		#if except_network_flag
		#	setting_expdf[!,:expvarall] = exp.((X_base[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
		#end

		# Update choice probabilities
		setting_expsum = combine(groupby(setting_expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
		if inattention_flag & (year > 2014)
			entrants = findall(households[!,:previous_plan_number] .== "NA");
			search_probabilities = 1 ./ (1 .+ exp.(Xsearch_base * beta_search));
			search_probabilities[entrants] .= 1; 
					
			setting_uninsured_probabilities = 1  ./ (1 .+ setting_expsum);
			setting_probabilities = search_probabilities[household_numbers] .* setting_expdf[!,:expvarall] ./ (1 .+ setting_expsum[household_numbers]) .+ 
				(1 .- search_probabilities[household_numbers]) .* data[uninsured_plans .== 0,:previous_choice];
		else
			setting_uninsured_probabilities = 1  ./ (1 .+ (setting_expsum).^(lambda));
			setting_probabilities = setting_expdf[!,:expvarall] .* (setting_expsum[household_numbers]).^(lambda-1) ./ (1 .+ (setting_expsum[household_numbers]).^(lambda));
		end
		setting_probability_weights = setting_probabilities .* hdata[:,findall(hdata_fields .== "weight")[1]];
		
		setting_all_probabilities = zeros(size(X_base)[1]);
		setting_all_probabilities[uninsured_plans .== 0] = setting_probabilities;
		setting_all_probabilities[uninsured_plans .== 1] = setting_uninsured_probabilities;
		setting_all_probability_weights = setting_all_probabilities .* data[!,:weight];
		
		setting_im_weight_object = DataFrame(household_number = household_numbers_data,
			bronze_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== .6),
			silver_weights = setting_all_probability_weights .* ((X_base[:,findall(all_covariates .== "AV")[1]] .== .7) .|
														   (X_base[:,findall(all_covariates .== "AV")[1]] .== .73) .|
														   (X_base[:,findall(all_covariates .== "AV")[1]] .== .87) .|
														   (X_base[:,findall(all_covariates .== "AV")[1]] .== .94)),
			gold_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== .8),
			platinum_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== .9),
			uninsured_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "AV")[1]] .== 0),
			anthem_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "Anthem")[1]] .== 1),
			bs_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "Blue_Shield")[1]] .== 1),
			hn_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "Health_Net")[1]] .== 1),
			kaiser_weights = setting_all_probability_weights .* (X_base[:,findall(all_covariates .== "Kaiser")[1]] .== 1),
			small_weights = setting_all_probability_weights .* (1 .- (X_base[:,findall(all_covariates .== "Anthem")[1]] .+
				X_base[:,findall(all_covariates .== "Blue_Shield")[1]] .+
				X_base[:,findall(all_covariates .== "Health_Net")[1]] .+ 
				X_base[:,findall(all_covariates .== "Kaiser")[1]])) .*  
				(X_base[:,findall(all_covariates .== "AV")[1]] .> 0));
		
		setting_im_household_probability_weights = 
			DataFrame(Bronze = combine(groupby(setting_im_weight_object,:household_number), :bronze_weights => sum)[!,:bronze_weights_sum],
					  Silver = combine(groupby(setting_im_weight_object,:household_number), :silver_weights => sum)[!,:silver_weights_sum],
					  Gold = combine(groupby(setting_im_weight_object,:household_number), :gold_weights => sum)[!,:gold_weights_sum],
					  Platinum = combine(groupby(setting_im_weight_object,:household_number), :platinum_weights => sum)[!,:platinum_weights_sum],
					  Uninsured = combine(groupby(setting_im_weight_object,:household_number), :uninsured_weights => sum)[!,:uninsured_weights_sum],
					  Anthem = combine(groupby(setting_im_weight_object,:household_number), :anthem_weights => sum)[!,:anthem_weights_sum],
					  Blue_Shield = combine(groupby(setting_im_weight_object,:household_number), :bs_weights => sum)[!,:bs_weights_sum],
					  Health_Net = combine(groupby(setting_im_weight_object,:household_number), :hn_weights => sum)[!,:hn_weights_sum],
					  Kaiser = combine(groupby(setting_im_weight_object,:household_number), :kaiser_weights => sum)[!,:kaiser_weights_sum],
					  Small_Insurer = combine(groupby(setting_im_weight_object,:household_number), :small_weights => sum)[!,:small_weights_sum]);
			

	# Get Scenario probabilities	
		
		# Initialize key variables
		subsidy = copy(voucher);
		new_unsubs_premiums = base_premium[household_pt_indices] .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		X_demand = copy(X);
		if inattention_flag
			Xsearch_df = copy(Xsearch);
		end
		
		# Create data frame to hold objects we will group by household
		silverdf = DataFrame(household_number = household_numbers[silver_indicators],premium = new_unsubs_premiums[silver_indicators],silver_indices = silver_indicators);
		
		# Update Subsidy
		
			# Determine SLC plans by household
			num_household_silver_plans = combine(groupby(silverdf,:household_number), :household_number => length)[!,:household_number_length];
			pop!(num_household_silver_plans)
			num_household_silver_plans = cat(dims=1,0,num_household_silver_plans);
			starting_indices = cumsum(num_household_silver_plans); # starting position for each household
			household_slc_indices = combine(groupby(silverdf,:household_number), :premium => determine_SLC_plan)[!,:premium_determine_SLC_plan];
			slc_indices = silverdf[starting_indices .+ household_slc_indices,:silver_indices]; # SLC indices in hdata
			slc_plans = zeros(size(hdata)[1]);
			slc_plans[slc_indices] = (households[!,:subsidized_members] .> 0); # SLC index gets 1 if house is subsidized, 0 otherwise
		
			# Determine new subsidy (unless voucher)
			if !voucher_flag
				subsidy = max.(0,new_unsubs_premiums[slc_indices] - households[!,:SLC_contribution]) .* 
					(households[!,:subsidized_members] .> 0); # each household's subsidy (0 if unsubsidized);
				subsidy = subsidy[household_numbers];
			end
			
			net_subsidies = min.(subsidy,new_unsubs_premiums);
			
		# Update Premium Variables
		X_demand[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]] = 
				(max.(0,new_unsubs_premiums .- subsidy) .- 
				(1 - mand_repeal_flag) .* hdata[:,findall(hdata_fields .== "penalty")[1]]) ./ hdata[:,findall(hdata_fields .== "members")[1]];			
		X_demand[uninsured_plans .== 0,premium_indices] = 
				X_demand[uninsured_plans .== 0,premium_indices] .* X_demand[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]];
		if random_parameters
			X_demand[:,random_parameter_indices] = X_demand[:,random_parameter_indices] .* V[:,:,1]; 
		end
		
		if inattention_flag & (year > 2014) & search_flag
		
			# Growth rates
			growth_rates = new_unsubs_premiums ./ forward_simulated_prev_premiums[uninsured_plans .== 0];
			fsim_prev_choices = forward_simulated_prev_choices[uninsured_plans .== 0];
			growthdf = DataFrame(household_number = household_numbers[offered_indices],growth_rate = growth_rates[offered_indices],prev_choice=fsim_prev_choices[offered_indices]);
			median_growth_rates = combine(groupby(growthdf,:household_number), :growth_rate => median);
			
			household_median_growth_rates = zeros(size(households)[1]);
			household_median_growth_rates[median_growth_rates[!,:household_number]] = median_growth_rates[!,:growth_rate_median];
			
			growthdf_for_choices = growthdf[findall(growthdf[!,:prev_choice] .== 1),:];
			growthdf_for_choices[!,:median_growth_rate] = household_median_growth_rates[growthdf_for_choices[!,:household_number]];
			
			premium_shock_households = growthdf_for_choices[findall(growthdf_for_choices[!,:growth_rate] .> growthdf_for_choices[!,:median_growth_rate]),:household_number];
			Xsearch_df[:,findall(search_covariates .== "premium_shock")[1]] .= 0;
			Xsearch_df[premium_shock_households,findall(search_covariates .== "premium_shock")[1]] .= 1;
		
		end
		
		expdf = DataFrame(household_number = household_numbers,
					expvarall = exp.((X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator .* no_inertia_flag)))/lambda));
		if except_network_flag
			expdf[!,:expvarall] = exp.((X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator_net .* except_network_flag)))/lambda);
		end
		
		# Update choice probabilities
		expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
		
		if inattention_flag & search_flag & (year > 2014)
			entrants = findall(households[!,:previous_plan_number] .== "NA");
			search_probabilities = 1 ./ (1 .+ exp.(Xsearch_df * beta_search));
			search_probabilities[entrants] .= 1; 
				
			uninsured_probabilities = 1  ./ (1 .+ expsum);
			probabilities = search_probabilities[household_numbers] .* expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers]) .+ 
				(1 .- search_probabilities[household_numbers]) .* data[uninsured_plans .== 0,:previous_choice];
		else
			uninsured_probabilities = 1  ./ (1 .+ (expsum).^(lambda));
			probabilities = expdf[!,:expvarall] .* (expsum[household_numbers]).^(lambda-1) ./ (1 .+ (expsum[household_numbers]).^(lambda));
		end
		
		all_probabilities = zeros(size(X_demand)[1]);
		all_probabilities[uninsured_plans .== 0] = probabilities;
		all_probabilities[uninsured_plans .== 1] = uninsured_probabilities;
		all_probability_weights = all_probabilities .* data[!,:weight];
		
		if no_uninsured_flag
			probabilities = probabilities ./ (1 .- uninsured_probabilities[household_numbers]);
			uninsured_probabilities .= 0;
		end
		
		# Compute consumer welfare
		
			# To compute consumer's surplus, we need utilities without penalty
			X_demand[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]] = 
				(max.(0,new_unsubs_premiums .- subsidy)) ./ hdata[:,findall(hdata_fields .== "members")[1]];					
			X_demand[uninsured_plans .== 0,premium_indices] = 
				X[uninsured_plans .== 0,premium_indices] .* X[uninsured_plans .== 0,findall(all_covariates .== "Premium")[1]];
					
			cwdf = DataFrame(household_number = household_numbers,
				alpha = hdata[:,findall(hdata_fields .== "alpha")[1]],
				U_incorrect = X_demand[uninsured_plans .== 0,:] * (min_beta .* (1 .- previous_choice_indicator .* no_inertia_flag)),
				inertia_parameter = (1 - no_inertia_flag) .* (X_demand[uninsured_plans .== 0,previous_choice_indices] * min_beta[previous_choice_indices]));
			cwdf[!,:expvarall] = exp.(cwdf[!,:U_incorrect]./lambda);
			if !inattention_flag
				cwdf[!,:inertia_penalty] = -X_demand[uninsured_plans .== 0,findall(all_covariates .== "previous_choice")[1]] .* probabilities .* 
					(cwdf[!,:inertia_parameter] ./ cwdf[!,:alpha]);
			end
			cwdf[!,:outside_utility] = cwdf[!,:alpha] .* hdata[:,findall(hdata_fields .== "penalty")[1]];
			
			alpha_household = combine(groupby(cwdf,:household_number), :alpha => minimum)[!,:alpha_minimum];
			U_outside = combine(groupby(cwdf,:household_number), :outside_utility => minimum)[!,:outside_utility_minimum];
			if inattention_flag
				inertia_penalty = 0
			else
				inertia_penalty = combine(groupby(cwdf,:household_number), :inertia_penalty => sum)[!,:inertia_penalty_sum];
			end
			expsum = combine(groupby(cwdf,:household_number), :expvarall => sum)[!,:expvarall_sum];	
			
			cw = -(1 ./ alpha_household) .* log.(expsum.^(lambda) + exp.(U_outside)) .- inertia_penalty;
			cw0 = -(1 ./ alpha_household) .* log.(expsum.^(lambda) + exp.(U_outside))
			cw75 = .75 * cw + .25 * cw0;
			cw50 = .50 * cw + .50 * cw0;
			cw25 = .25 * cw + .75 * cw0;
			
		# Intensive margin probabilities
		expdf[!,:probabilities] = probabilities;
		expdf[!,:bronze_probabilities] = probabilities .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .6)
		expdf[!,:silver_probabilities] = probabilities .* ((X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .7) .|
														   (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .73) .|
														   (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .87) .|
														   (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .94))
		expdf[!,:gold_probabilities] = probabilities .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .8)
		expdf[!,:platinum_probabilities] = probabilities .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .9)
		
		expdf[!,:aca_probabilities] = aca_probabilities;
		expdf[!,:bronze_aca_probabilities] = aca_probabilities .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .6)
		expdf[!,:silver_aca_probabilities] = aca_probabilities .* ((X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .7) .|
														   (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .73) .|
														   (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .87) .|
														   (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .94))
		expdf[!,:gold_aca_probabilities] = aca_probabilities .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .8)
		expdf[!,:platinum_aca_probabilities] = aca_probabilities .* (X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .9)
		
		probs_household = combine(groupby(expdf,:household_number), :probabilities => sum)[!,:probabilities_sum];
		aca_probs_household = combine(groupby(expdf,:household_number), :aca_probabilities => sum)[!,:aca_probabilities_sum];
		intensive_probabilities = probabilities .* (combine(groupby(expdf,:household_number), :aca_probabilities => sum)[!,:aca_probabilities_sum] ./ 
			combine(groupby(expdf,:household_number), :probabilities => sum)[!,:probabilities_sum])[household_numbers];

		X_bronze_indices = findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .6);
		X_silver_indices = cat(dims=1,findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .7),
					findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .73),
					findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .87),
					findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .94));
		X_brsil_indices = cat(dims=1,X_bronze_indices,X_silver_indices);
		X_gold_indices = findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .8);
		X_platinum_indices = findall(X_demand[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]] .== .9);

		intensive_metal_probabilities = zeros(size(hdata)[1]);
		intensive_metal_probabilities[X_bronze_indices] = probabilities[X_bronze_indices] .* 
			((combine(groupby(expdf,:household_number), :bronze_aca_probabilities => sum)[!,:bronze_aca_probabilities_sum] ./ 
			combine(groupby(expdf,:household_number), :bronze_probabilities => sum)[!,:bronze_probabilities_sum])[household_numbers])[X_bronze_indices];
		intensive_metal_probabilities[X_silver_indices] = probabilities[X_silver_indices] .* 
			((combine(groupby(expdf,:household_number), :silver_aca_probabilities => sum)[!,:silver_aca_probabilities_sum] ./ 
			combine(groupby(expdf,:household_number), :silver_probabilities => sum)[!,:silver_probabilities_sum])[household_numbers])[X_silver_indices];
		intensive_metal_probabilities[X_gold_indices] = probabilities[X_gold_indices] .* 
			((combine(groupby(expdf,:household_number), :gold_aca_probabilities => sum)[!,:gold_aca_probabilities_sum] ./ 
			combine(groupby(expdf,:household_number), :gold_probabilities => sum)[!,:gold_probabilities_sum])[household_numbers])[X_gold_indices];
		intensive_metal_probabilities[X_platinum_indices] = probabilities[X_platinum_indices] .* 
			((combine(groupby(expdf,:household_number), :platinum_aca_probabilities => sum)[!,:platinum_aca_probabilities_sum] ./ 
			combine(groupby(expdf,:household_number), :platinum_probabilities => sum)[!,:platinum_probabilities_sum])[household_numbers])[X_platinum_indices];
		
		# Prepare CSR variables - factors to apply
			# 73% AV plan - government covers 3%, no moral hazard assumed in CMS paper
			# 87% AV plan - government covers 17%, 12% moral hazard assumed in CMS paper
			# 94% AV plan - government covers 24%, 12% moral hazard assumed in CMS paper
		silver_avs = X[uninsured_plans .== 0,findall(all_covariates .== "AV")[1]];
		silver_avs[silver_avs .== .6] .= 0;
		silver_avs[silver_avs .== .8] .= 0;
		silver_avs[silver_avs .== .9] .= 0;
		government_share = zeros(length(silver_avs));
		government_share[findall(silver_avs .== 0.73)] .= 0.03;
		government_share[findall(silver_avs .== 0.87)] .= 0.17 * 1.12;
		government_share[findall(silver_avs .== 0.94)] .= 0.24 * 1.12;
		
		# Calculate weights
		weights = hdata[:,findall(hdata_fields .== "weight")[1]];
		probability_weights = probabilities .* weights;
		intensive_probability_weights = intensive_probabilities .* weights;
		intensive_metal_probability_weights = intensive_metal_probabilities .* weights;
		pricing_weights = probabilities .* hdata[:,findall(hdata_fields .== "age_factor")[1]];
		aca_pricing_weights = aca_probabilities .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		rf_pricing_weights = probabilities .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		intensive_pricing_weights = intensive_probabilities .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		intensive_metal_pricing_weights = intensive_metal_probabilities .* hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		
		if all_rs_variables
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names, plan_pt_name = data_pt_names, weights=weights,
				probability_weights = probability_weights,intensive_probability_weights = intensive_probability_weights,intensive_metal_probability_weights = intensive_metal_probability_weights,
				pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,aca_pricing_weights = aca_pricing_weights,
					intensive_pricing_weights = intensive_pricing_weights, intensive_metal_pricing_weights = intensive_metal_pricing_weights,
				net_subsidies = net_subsidies .* probabilities, government_share_weights = government_share .* weights,
				rs_share_weights_0to34 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]],
				rs_share_weights_35to54 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_35to54")[1]],
				rs_share_weights_250to400 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_250to400")[1]],
				rs_share_weights_gt400 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_gt400")[1]],
				rs_share_weights_male = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]],
				rs_share_weights_family = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]],
				rs_share_weights_Asian = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Asian")[1]],
				rs_share_weights_Black = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Black")[1]],
				rs_share_weights_Hispanic = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]],
				rs_share_weights_Other = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Other_Race")[1]]);
		elseif int_rs_variables
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names, plan_pt_name = data_pt_names, weights=weights,
				probability_weights = probability_weights,intensive_probability_weights = intensive_probability_weights,intensive_metal_probability_weights = intensive_metal_probability_weights,
				pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,aca_pricing_weights = aca_pricing_weights,
					intensive_pricing_weights = intensive_pricing_weights, intensive_metal_pricing_weights = intensive_metal_pricing_weights,
				net_subsidies = net_subsidies .* probabilities, government_share_weights = government_share .* weights,
				rs_share_weights_0to34 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]],
				rs_share_weights_male = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]],
				rs_share_weights_family = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]],
				rs_share_weights_minority = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_minority")[1]]);
		elseif old_flag
			weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names, plan_pt_name = data_pt_names, weights=weights,
				probability_weights = probability_weights,intensive_probability_weights = intensive_probability_weights,intensive_metal_probability_weights = intensive_metal_probability_weights,
				pricing_weights = pricing_weights,rf_pricing_weights = rf_pricing_weights,aca_pricing_weights = aca_pricing_weights,
					intensive_pricing_weights = intensive_pricing_weights, intensive_metal_pricing_weights = intensive_metal_pricing_weights,
				net_subsidies = net_subsidies .* probabilities, government_share_weights = government_share .* weights,
				rs_share_weights_18to54 = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_18to54")[1]],
				rs_share_weights_Hispanic = probability_weights .* rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]]);
		end
		
		
		im_weight_object = DataFrame(household_number = household_numbers_data,
			bronze_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "AV")[1]] .== .6),
			silver_weights = all_probability_weights .* ((X_demand[:,findall(all_covariates .== "AV")[1]] .== .7) .|
														   (X_demand[:,findall(all_covariates .== "AV")[1]] .== .73) .|
														   (X_demand[:,findall(all_covariates .== "AV")[1]] .== .87) .|
														   (X_demand[:,findall(all_covariates .== "AV")[1]] .== .94)),
			gold_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "AV")[1]] .== .8),
			platinum_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "AV")[1]] .== .9),
			uninsured_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "AV")[1]] .== 0),
			anthem_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "Anthem")[1]] .== 1),
			bs_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "Blue_Shield")[1]] .== 1),
			hn_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "Health_Net")[1]] .== 1),
			kaiser_weights = all_probability_weights .* (X_demand[:,findall(all_covariates .== "Kaiser")[1]] .== 1),
			small_weights = all_probability_weights .* (1 .- (X_demand[:,findall(all_covariates .== "Anthem")[1]] .+
				X_demand[:,findall(all_covariates .== "Blue_Shield")[1]] .+
				X_demand[:,findall(all_covariates .== "Health_Net")[1]] .+ 
				X_demand[:,findall(all_covariates .== "Kaiser")[1]])) .*  
				(X_demand[:,findall(all_covariates .== "AV")[1]] .> 0));
		
		im_household_probability_weights = 
			DataFrame(Bronze = combine(groupby(im_weight_object,:household_number), :bronze_weights => sum)[!,:bronze_weights_sum],
					  Silver = combine(groupby(im_weight_object,:household_number), :silver_weights => sum)[!,:silver_weights_sum],
					  Gold = combine(groupby(im_weight_object,:household_number), :gold_weights => sum)[!,:gold_weights_sum],
					  Platinum = combine(groupby(im_weight_object,:household_number), :platinum_weights => sum)[!,:platinum_weights_sum],
					  Uninsured = combine(groupby(im_weight_object,:household_number), :uninsured_weights => sum)[!,:uninsured_weights_sum],
					  Anthem = combine(groupby(im_weight_object,:household_number), :anthem_weights => sum)[!,:anthem_weights_sum],
					  Blue_Shield = combine(groupby(im_weight_object,:household_number), :bs_weights => sum)[!,:bs_weights_sum],
					  Health_Net = combine(groupby(im_weight_object,:household_number), :hn_weights => sum)[!,:hn_weights_sum],
					  Kaiser = combine(groupby(im_weight_object,:household_number), :kaiser_weights => sum)[!,:kaiser_weights_sum],
					  Small_Insurer = combine(groupby(im_weight_object,:household_number), :small_weights => sum)[!,:small_weights_sum]);
			
		# Update plan-level variables
		plan_pricing_weights = sort!(combine(groupby(weight_object,:plan_name), :pricing_weights => sum))[!,:pricing_weights_sum];
		plan_pt_demands = sort!(combine(groupby(weight_object,:plan_pt_name), :probability_weights => sum))[!,:probability_weights_sum];
		plan_demands = sort!(combine(groupby(weight_object,:plan_name), :probability_weights => sum))[!,:probability_weights_sum];
		intensive_plan_demands = sort!(combine(groupby(weight_object,:plan_name), :intensive_probability_weights => sum))[!,:intensive_probability_weights_sum];
		intensive_metal_plan_demands = sort!(combine(groupby(weight_object,:plan_name), :intensive_metal_probability_weights => sum))[!,:intensive_metal_probability_weights_sum];
		plan_total_premiums = base_premium[year_plan_indices] .* 
			sort!(combine(groupby(weight_object,:plan_name), :rf_pricing_weights => sum))[!,:rf_pricing_weights_sum];
		aca_plan_total_premiums = base_premium[year_plan_indices] .* 
			sort!(combine(groupby(weight_object,:plan_name), :aca_pricing_weights => sum))[!,:aca_pricing_weights_sum];
		intensive_plan_total_premiums = base_premium[year_plan_indices] .* 
			sort!(combine(groupby(weight_object,:plan_name), :intensive_pricing_weights => sum))[!,:intensive_pricing_weights_sum];
		intensive_metal_plan_total_premiums = base_premium[year_plan_indices] .* 
			sort!(combine(groupby(weight_object,:plan_name), :intensive_metal_pricing_weights => sum))[!,:intensive_metal_pricing_weights_sum];
		plan_subsidies = sort!(combine(groupby(weight_object,:plan_name), :net_subsidies => sum))[!,:net_subsidies_sum];
		plan_government_share = sort!(combine(groupby(weight_object,:plan_name), :government_share_weights => sum))[!,:government_share_weights_sum] ./
			sort!(combine(groupby(weight_object,:plan_name), :weights => sum))[!,:weights_sum];
		
		# Update pmt variables
		X_RS_new = copy(X_RS);
		if all_rs_variables
			X_RS_new[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_0to34 => sum),[:plan_name])[!,:rs_share_weights_0to34_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_35to54")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_35to54 => sum),[:plan_name])[!,:rs_share_weights_35to54_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_250to400")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_250to400 => sum),[:plan_name])[!,:rs_share_weights_250to400_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_gt400")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_gt400 => sum),[:plan_name])[!,:rs_share_weights_gt400_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_male => sum),[:plan_name])[!,:rs_share_weights_male_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_family => sum),[:plan_name])[!,:rs_share_weights_family_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Asian")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Asian => sum),[:plan_name])[!,:rs_share_weights_Asian_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Black")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Black => sum),[:plan_name])[!,:rs_share_weights_Black_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Hispanic => sum),[:plan_name])[!,:rs_share_weights_Hispanic_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Other_Race")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Other => sum),[:plan_name])[!,:rs_share_weights_Other_sum] ./ plan_demands;
		elseif int_rs_variables
			X_RS_new[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_0to34 => sum),[:plan_name])[!,:rs_share_weights_0to34_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_male => sum),[:plan_name])[!,:rs_share_weights_male_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_family => sum),[:plan_name])[!,:rs_share_weights_family_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_minority")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_minority => sum),[:plan_name])[!,:rs_share_weights_minority_sum] ./ plan_demands;
		elseif old_flag
			X_RS_new[:,findall(risk_score_covariates .== "share_18to54")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_18to54 => sum),[:plan_name])[!,:rs_share_weights_18to54_sum] ./ plan_demands;
			X_RS_new[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :rs_share_weights_Hispanic => sum),[:plan_name])[!,:rs_share_weights_Hispanic_sum] ./ plan_demands;
		end
		
		plan_risk_scores = exp.(X_RS_new * theta_RS);
		household_risk_scores = plan_risk_scores[household_pmt_indices];
		
		# Calculate RA Transfers
		plan_transfers = zeros(length(plan_pmt_names));
		if !nora_flag
			plan_rs_shares = plan_demands .* plan_risk_scores ./ 
				sum(plan_demands .* plan_risk_scores);
			plan_ra_shares = plan_demands .* pmt[:,findall(pmt_fields .== "ra_factor")[1]] ./ 
				sum(plan_demands .* pmt[:,findall(pmt_fields .== "ra_factor")[1]]);
			plan_transfers = (plan_rs_shares .- plan_ra_shares) .* sum(plan_total_premiums);
		end
		
		# Calculate Average Claims
		X_CLMS_new = X_CLMS[1:size(X_CLMS)[1],:];
		X_CLMS_new[:,findall(claims_moments_covariates .== "log_risk_score")[1]] = X_RS_new * theta_RS;
		average_claims = exp.(X_CLMS_new * theta_CLMS);	
		average_pt_claims = (transpose(foc_moment_indicator) * (average_claims .* plan_demands)) ./
					(transpose(foc_moment_indicator) * plan_demands);
		
		# Calculate Profit and FOC Residuals
		foc_resids = update_foc_values(base_premium,pmt,hdata,X_RS,X_CLMS,foc_residuals,nora_flag,mand_repeal_flag,voucher_flag,pc_flag,no_inertia_flag,no_uninsured_flag,ef_approach_flag);
		foc_resids_by_metal = calculate_foc_resids_by_metal(plan_pt_demands,foc_resids);
		foc_resids_by_firm = calculate_foc_resids_by_firm(plan_pt_demands,foc_resids);
		plan_profits = compute_plan_profit(plan_total_premiums,average_claims,plan_transfers,plan_demands,false);
		output[findall(output_fields .== "Average FOC Residuals")[1]] = sum(plan_pt_demands .* foc_resids)/sum(plan_pt_demands);
		output[findall(output_fields .== "Max FOC Residuals")[1]] = maximum(abs.(foc_resids));
		
		# Uncompensated Care Calculation:
			# I got the 1942 number from an Urban/Kaiser report (https://www.kff.org/uninsured/report/uncompensated-care-for-the-uninsured-in-2013-a-detailed-examination/)
			# I got the 4.3% inflation factor from CMS/NHEA data (https://www.cms.gov/research-statistics-data-and-systems/statistics-trends-and-reports/nationalhealthexpenddata/nhe-fact-sheet.html) 
		avg_uncomp_care = (1942/12) * 1.043;
		
		# Welfare fields
		output[findall(output_fields .== "Profit")[1]] = sum(plan_profits .* plan_pt_demands)/total_weight;
		output[findall(output_fields .== "Consumer Surplus")[1]] = sum(cw .* households[!,:weight]) / total_weight;
		output[findall(output_fields .== "Premium Subsidies")[1]] = sum(net_subsidies .* probabilities) / total_weight;
		output[findall(output_fields .== "Cost Sharing Subsidies")[1]] = 
			sum(plan_government_share .* average_claims .* plan_demands) / total_weight;
		output[findall(output_fields .== "Penalties")[1]] = (1 - mand_repeal_flag) * 
			sum(uninsured_probabilities .* (households[!,:penalty] ./ 12) .* households[!,:weight]) / total_weight;	
	 	if !no_uninsured_flag
			output[findall(output_fields .== "Uncompensated Care")[1]] =
				sum(uninsured_probabilities .* avg_uncomp_care .* households[!,:weight]) / total_weight;	
		end
		output[findall(output_fields .== "Social Welfare")[1]] =
			output[findall(output_fields .== "Profit")[1]] + 
			output[findall(output_fields .== "Consumer Surplus")[1]] + 
			output[findall(output_fields .== "Penalties")[1]] - 
			output[findall(output_fields .== "Premium Subsidies")[1]] - 
			output[findall(output_fields .== "Cost Sharing Subsidies")[1]] - 
			output[findall(output_fields .== "Uncompensated Care")[1]];
		
		
		output[findall(output_fields .== "CWlt250")[1]] = sum(cw[findall(households[!,:FPL] .<= 2.5)] .* 
				households[findall(households[!,:FPL] .<= 2.5),:weight]) / 
			sum(households[findall(households[!,:FPL] .<= 2.5),:weight]);
		output[findall(output_fields .== "CW250to400")[1]] = sum(cw[findall((households[!,:FPL] .> 2.50) .& (households[!,:FPL] .<= 4))] .* 
				households[findall((households[!,:FPL] .> 2.50) .& (households[!,:FPL] .<= 4)),:weight]) / 
			sum(households[findall((households[!,:FPL] .> 2.50) .& (households[!,:FPL] .<= 4)),:weight]);
		output[findall(output_fields .== "CWgt400")[1]] = sum(cw[findall(households[!,:FPL] .> 4)] .*
				households[findall(households[!,:FPL] .> 4),:weight]) / 
			sum(households[findall(households[!,:FPL] .> 4),:weight]);
		output[findall(output_fields .== "CW0to17")[1]] = 
			sum(cw .* households[!,:perc_0to17] .* households[!,:weight]) / 
			sum(households[!,:perc_0to17] .* households[!,:weight]);
		output[findall(output_fields .== "CW18to34")[1]] = 
			sum(cw .* (households[!,:perc_18to25] .+ households[!,:perc_26to34]) .* households[!,:weight]) / 
			sum((households[!,:perc_18to25] .+ households[!,:perc_26to34]) .* households[!,:weight]);
		output[findall(output_fields .== "CW35to54")[1]] = 
			sum(cw .* (households[!,:perc_35to44] .+ households[!,:perc_45to54]) .* households[!,:weight]) / 
			sum((households[!,:perc_35to44] .+ households[!,:perc_45to54]) .* households[!,:weight]);
		output[findall(output_fields .== "CW55plus")[1]] = 
			sum(cw .* (households[!,:perc_55to64] .+ households[!,:perc_65plus]) .* households[!,:weight]) / 
			sum((households[!,:perc_55to64] .+ households[!,:perc_65plus]) .* households[!,:weight]);
		output[findall(output_fields .== "CWAsian")[1]] = 
			sum(cw .* households[!,:perc_asian] .* households[!,:weight]) / 
			sum(households[!,:perc_asian] .* households[!,:weight]);
		output[findall(output_fields .== "CWBlack")[1]] = 
			sum(cw .* households[!,:perc_black] .* households[!,:weight]) / 
			sum(households[!,:perc_black] .* households[!,:weight]);
		output[findall(output_fields .== "CWHispanic")[1]] = 
			sum(cw .* households[!,:perc_hispanic] .* households[!,:weight]) / 
			sum(households[!,:perc_hispanic] .* households[!,:weight]);
		output[findall(output_fields .== "CWWhite")[1]] = 
			sum(cw .* households[!,:perc_white] .* households[!,:weight]) / 
			sum(households[!,:perc_white] .* households[!,:weight]);
		output[findall(output_fields .== "CWOther")[1]] = 
			sum(cw .* households[!,:perc_other] .* households[!,:weight]) / 
			sum(households[!,:perc_other] .* households[!,:weight]);
		output[findall(output_fields .== "CWMale")[1]] = 
			sum(cw .* households[!,:perc_male] .* households[!,:weight]) / 
			sum(households[!,:perc_male] .* households[!,:weight]);
		output[findall(output_fields .== "CWFemale")[1]] = 
			sum(cw .* (1 .- households[!,:perc_male]) .* households[!,:weight]) / 
			sum((1 .- households[!,:perc_male]) .* households[!,:weight]);
		output[findall(output_fields .== "CWSingle")[1]] = 
			sum(cw[findall(households[!,:household_size] .== 1)] .* households[findall(households[!,:household_size] .== 1),:weight]) / 
			sum(households[findall(households[!,:household_size] .== 1),:weight]);
		output[findall(output_fields .== "CWFamily")[1]] = 
			sum(cw[findall(households[!,:household_size] .> 1)] .* households[findall(households[!,:household_size] .> 1),:weight]) / 
			sum(households[findall(households[!,:household_size] .> 1),:weight]);
		
		
		
		output[findall(output_fields .== "Consumer Surplus - Kappa 75%")[1]] = sum(cw75 .* households[!,:weight]) / total_weight;
		output[findall(output_fields .== "Consumer Surplus - Kappa 50%")[1]] = sum(cw50 .* households[!,:weight]) / total_weight;
		output[findall(output_fields .== "Consumer Surplus - Kappa 25%")[1]] = sum(cw25 .* households[!,:weight]) / total_weight;
		output[findall(output_fields .== "Consumer Surplus - Kappa 0%")[1]] = sum(cw0 .* households[!,:weight]) / total_weight;
		
		# Coverage Fields
		output[findall(in(["Bronze","Silver","Gold","Platinum","Total Coverage"]),output_fields)] = 
			tabulate_coverage_by_metal(plan_demands);
		output[findall(in(["Bronze","Silver","Gold","Platinum"]),output_fields)] += get_dominated_enrollment(probability_weights);
		output[findall(in(["HMO","PPO"]),output_fields)] = tabulate_coverage_by_network_type(plan_demands);
		output[findall(in(insurers),output_fields)] = tabulate_coverage_by_insurer(plan_demands);
		output[findall(output_fields .== "% Enrolled")[1]] = output[findall(output_fields .== "Total Coverage")[1]] / total_pop;
		
		# Switching Rates
		output[findall(output_fields .== "% Switching Plan (Base)")[1]] = 
			#(sum(abs.(base_probability_weights .- probability_weights)) / 2)/sum(base_probability_weights);
			(sum(abs.(base_all_probability_weights .- all_probability_weights)) / 2)/sum(households[!,:weight]);
		output[findall(output_fields .== "% Switching Plan (Setting)")[1]] = 
			#(sum(abs.(setting_probability_weights .- probability_weights)) / 2)/sum(setting_probability_weights);
			(sum(abs.(setting_all_probability_weights .- all_probability_weights)) / 2)/sum(households[!,:weight]);
		output[findall(output_fields .== "% Switching Metal (Base)")[1]] = 
			(sum(abs.(base_im_household_probability_weights[!,:Bronze] .- im_household_probability_weights[!,:Bronze])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Silver] .- im_household_probability_weights[!,:Silver])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Gold] .- im_household_probability_weights[!,:Gold])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Platinum] .- im_household_probability_weights[!,:Platinum])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Uninsured] .- im_household_probability_weights[!,:Uninsured])))/
			(2 .* sum(households[!,:weight]));
		output[findall(output_fields .== "% Switching Metal (Setting)")[1]] = 
			(sum(abs.(setting_im_household_probability_weights[!,:Bronze] .- im_household_probability_weights[!,:Bronze])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Silver] .- im_household_probability_weights[!,:Silver])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Gold] .- im_household_probability_weights[!,:Gold])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Platinum] .- im_household_probability_weights[!,:Platinum])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Uninsured] .- im_household_probability_weights[!,:Uninsured])))/
			(2 .* sum(households[!,:weight]));
		output[findall(output_fields .== "% Switching Insurer (Base)")[1]] = 
			(sum(abs.(base_im_household_probability_weights[!,:Anthem] .- im_household_probability_weights[!,:Anthem])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Blue_Shield] .- im_household_probability_weights[!,:Blue_Shield])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Health_Net] .- im_household_probability_weights[!,:Health_Net])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Kaiser] .- im_household_probability_weights[!,:Kaiser])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Small_Insurer] .- im_household_probability_weights[!,:Small_Insurer])) .+ 
			sum(abs.(base_im_household_probability_weights[!,:Uninsured] .- im_household_probability_weights[!,:Uninsured])))/
			(2 .* sum(households[!,:weight]));
		output[findall(output_fields .== "% Switching Insurer (Setting)")[1]] = 
			(sum(abs.(setting_im_household_probability_weights[!,:Anthem] .- im_household_probability_weights[!,:Anthem])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Blue_Shield] .- im_household_probability_weights[!,:Blue_Shield])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Health_Net] .- im_household_probability_weights[!,:Health_Net])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Kaiser] .- im_household_probability_weights[!,:Kaiser])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Small_Insurer] .- im_household_probability_weights[!,:Small_Insurer])) .+ 
			sum(abs.(setting_im_household_probability_weights[!,:Uninsured] .- im_household_probability_weights[!,:Uninsured])))/
			(2 .* sum(households[!,:weight]));
		
		# Premium Fields
		output[findall(in(["Bronze Premium";"Silver Premium";"Gold Premium";"Platinum Premium"]),output_fields)] =
			compute_avg_premium_by_metal(plan_total_premiums,plan_demands);
		output[findall(in(["Anthem Premium";"Blue Shield Premium";"Health Net Premium";"Kaiser Premium";"Other Insurer Premium"]),output_fields)] =
			compute_avg_premium_by_firm(plan_total_premiums,plan_demands);
		output[findall(in(["HMO Premium";"PPO Premium"]),output_fields)] =
			compute_avg_premium_by_nettype(plan_total_premiums,plan_demands);
		output[findall(output_fields .== "Average Premium")[1]] = sum(plan_total_premiums)/sum(plan_demands);
		output[findall(output_fields .== "Average Subsidized Premium")[1]] = sum(plan_total_premiums .- plan_subsidies)/sum(plan_demands);
		
		output[findall(in(["Bronze Premium (No Selection)";"Silver Premium (No Selection)";
			"Gold Premium (No Selection)";"Platinum Premium (No Selection)"]),output_fields)] =
			compute_avg_premium_by_metal(aca_plan_total_premiums,aca_plan_demands);
		output[findall(in(["Bronze Premium (No Metal Selection)";"Silver Premium (No Metal Selection)";
			"Gold Premium (No Metal Selection)";"Platinum Premium (No Metal Selection)"]),output_fields)] =
			compute_avg_premium_by_metal(intensive_metal_plan_total_premiums,intensive_metal_plan_demands);
		output[findall(in(["Bronze Premium (Intensive Selection)";"Silver Premium (Intensive Selection)";
			"Gold Premium (Intensive Selection)";"Platinum Premium (Intensive Selection)"]),output_fields)] =
			compute_avg_premium_by_metal(intensive_plan_total_premiums,intensive_plan_demands);
			
		output[findall(output_fields .== "Average Premium (No Selection)")[1]] = sum(aca_plan_total_premiums)/sum(aca_plan_demands);	
		output[findall(output_fields .== "Average Premium (No Metal Selection)")[1]] = sum(intensive_metal_plan_total_premiums)/sum(intensive_metal_plan_demands);	
		output[findall(output_fields .== "Average Premium (Intensive Selection)")[1]] = sum(intensive_plan_total_premiums)/sum(intensive_plan_demands);	
			
		# Claims Fields
		output[findall(in(["Bronze Avg Claims Resid";"Silver Avg Claims Resid";"Gold Avg Claims Resid";"Platinum Avg Claims Resid"]),output_fields)] =	
			compute_avg_claims_by_metal(average_claims,average_pt_claims,plan_demands) .+ foc_resids_by_metal;
		output[findall(in(["Bronze Avg Claims (No Metal Selection)";"Silver Avg Claims (No Metal Selection)";"Gold Avg Claims (No Metal Selection)";"Platinum Avg Claims (No Metal Selection)"]),output_fields)] =	
			compute_avg_claims_by_metal(average_claims,average_pt_claims,intensive_metal_plan_demands) .+ foc_resids_by_metal;
		output[findall(in(["Bronze Avg Claims (Intensive Selection)";"Silver Avg Claims (Intensive Selection)";"Gold Avg Claims (Intensive Selection)";"Platinum Avg Claims (Intensive Selection)"]),output_fields)] =	
			compute_avg_claims_by_metal(average_claims,average_pt_claims,intensive_plan_demands) .+ foc_resids_by_metal;
		output[findall(in(["Bronze Avg Claims (No Selection)";"Silver Avg Claims (No Selection)";"Gold Avg Claims (No Selection)";"Platinum Avg Claims (No Selection)"]),output_fields)] =
			compute_avg_claims_by_metal(aca_average_claims,average_pt_claims,plan_demands) .+ foc_resids_by_metal;		
		
		output[findall(in(["Anthem Avg Claims Resid";"Blue Shield Avg Claims Resid";"Health Net Avg Claims Resid";"Kaiser Avg Claims Resid";"Other Insurer Avg Claims Resid"]),output_fields)] =
			compute_avg_claims_by_insurer(average_claims,average_pt_claims,plan_demands) .+ foc_resids_by_firm;
		output[findall(in(["Anthem Avg Claims (No Selection)";"Blue Shield Avg Claims (No Selection)";"Health Net Avg Claims (No Selection)";"Kaiser Avg Claims (No Selection)";"Other Insurer Avg Claims (No Selection)"]),output_fields)] =
			compute_avg_claims_by_insurer(aca_average_claims,average_pt_claims,plan_demands) .+ foc_resids_by_firm;		
			
		output[findall(output_fields .== "Average Claims Resids")[1]] = 
			sum(plan_demands .* average_claims)/sum(plan_demands) .+
			output[findall(output_fields .== "Average FOC Residuals")[1]]; 
		output[findall(output_fields .== "Average Claims (No Metal Selection)")[1]] = 
			sum(intensive_metal_plan_demands .* average_claims)/sum(intensive_metal_plan_demands) .+
			output[findall(output_fields .== "Average FOC Residuals")[1]]; 
		output[findall(output_fields .== "Average Claims (Intensive Selection)")[1]] = 
			sum(intensive_plan_demands .* average_claims)/sum(intensive_plan_demands) .+
			output[findall(output_fields .== "Average FOC Residuals")[1]]; 
		output[findall(output_fields .== "Average Claims (No Selection)")[1]] = 
			sum(plan_demands .* aca_average_claims)/sum(plan_demands) .+
			output[findall(output_fields .== "Average FOC Residuals")[1]]; 
				
		output[findall(output_fields .== "Average Risk Score")[1]] = sum(plan_demands .* plan_risk_scores)/sum(plan_demands);
		output[findall(output_fields .== "Average AV")[1]] = sum(plan_demands .* plans_pmt[!,:AV])/sum(plan_demands);	
			
		if int_rs_variables
			uninsured_share_0to34 = sum(uninsured_probabilities .* (data[!,:weight] .* 
					(households[household_numbers_data,:perc_18to25] .+ households[household_numbers_data,:perc_26to34] .+
					households[household_numbers_data,:perc_35to44] .+ households[household_numbers_data,:perc_45to54]))[findall(uninsured_plans .== 1)])/
				sum(uninsured_probabilities .* data[findall(uninsured_plans .== 1),:weight]);
			uninsured_share_male = sum(uninsured_probabilities .* (data[!,:weight] .* households[household_numbers_data,:perc_male])[findall(uninsured_plans .== 1)])/
				sum(uninsured_probabilities .* data[findall(uninsured_plans .== 1),:weight]);	
			uninsured_share_family = sum(uninsured_probabilities .* (data[!,:weight] .* (households[household_numbers_data,:household_size] .> 1))[findall(uninsured_plans .== 1)])/
				sum(uninsured_probabilities .* data[findall(uninsured_plans .== 1),:weight]);	
			uninsured_share_minority = sum(uninsured_probabilities .* (data[!,:weight] .* (households[household_numbers_data,:perc_asian] .+ households[household_numbers_data,:perc_black] .+ 
					households[household_numbers_data,:perc_hispanic] .+ households[household_numbers_data,:perc_other]))[findall(uninsured_plans .== 1)])/
				sum(uninsured_probabilities .* data[findall(uninsured_plans .== 1),:weight]);	
			
			output[findall(output_fields .== "Uninsured Risk Score")[1]] =
				exp(sum(theta_RS .* cat(dims=1,1,0,0,0,uninsured_share_0to34,uninsured_share_male,uninsured_share_family,uninsured_share_minority)))
		end	
			
		#uninsured_share_18to54 = sum(uninsured_probabilities .* (data[!,:weight] .* 
		#		(households[household_numbers_data,:perc_18to25] .+ households[household_numbers_data,:perc_26to34] .+
		#		households[household_numbers_data,:perc_35to44] .+ households[household_numbers_data,:perc_45to54]))[findall(uninsured_plans .== 1)])/
		#	sum(uninsured_probabilities .* data[findall(uninsured_plans .== 1),:weight]);
		#uninsured_share_male = sum(uninsured_probabilities .* (data[!,:weight] .* households[household_numbers_data,:perc_male])[findall(uninsured_plans .== 1)])/
		#	sum(uninsured_probabilities .* data[findall(uninsured_plans .== 1),:weight]);	
		#output[findall(output_fields .== "Uninsured Risk Score")[1]] =
		#	exp(sum(theta_RS .* cat(dims=1,1,0,0,0,uninsured_share_18to54,uninsured_share_male)))
		
		# Scale uncompensated care by the % change in the uninsured risk score
		#if !no_uninsured_flag
		#	output[findall(output_fields .== "Uncompensated Care")[1]] = output[findall(output_fields .== "Uncompensated Care")[1]] *
		#		output[findall(output_fields .== "Uninsured Risk Score")[1]]/output[findall(output_fields .== "Uninsured Risk Score")[1],1];
		#end

	
	# Save output
	results = DataFrame(cat(dims=2,output_fields,output),:auto);
	cd(datadir)
	CSV.write(output_file,results);	
